function [Pd_patient_wise,FA_per_image,AUC,Pd,Pfa]=get_ROC_KDD(p, Y, patient_ID, FA_low, FA_high);
% [Pd_patient_wise,FA_per_image,AUC,Pd,Pfa]=get_ROC_KDD(p, Y, patient_ID, FP_low, FP_high);
% Compute: Malignant Patient detection rate <vs> False-alarmss-per-image 
%     FROC curves for KDD Cup 2008
% INPUTS: (the first 3 inputs have to be of the same size)                      
% p has predictions which are real numbers (ie not just in [0,1])
%    Larger p indicates higher confidence that label Y is positive.
% Y denotes positive class labels elements by +1.
%    All other values (eg 0, -1) in Y are considered -ve
% patient_ID has to be numbers, unique for each patient
% [FA_low,FA_high] is the region in which the area under FROC (AUC) has to be
% computed [default is [0.2,0.3]]
% OUTPUTS:
% Pd_patient_wise,FA_per_image: per-patient sensitivity and FA per image
%       (for the Official KDD Cup FROC)
% AUC is the area under the Official FROC (see above) in the FA range 
%       [FA_low, FA_high]
% Pd, Pfa: Probability of detection and probability of false alarm for
%       traditional ROC (not used for KDD Cup evaluation)
%------------------------------------------------------------
% Error checking:
%------------------------------------------------------------
if (rows(p)==1)
    p=p';
end;

if (rows(Y)==1)
    Y=Y';
end;

if (rows(patient_ID)==1)
    Patient_ID=Patient_ID';
end;

if (rows(p)~=rows(Y))
    error('p and Y are of different lengths');
end

if (rows(p)~=rows(patient_ID))
    error('p and patient_ID are of different lengths');
end

if ~exist('FA_low','var')
	FA_low=0.2;
end

if ~exist('FA_high','var')
	FA_high=0.3;
end
%------------------------------------------------------------
% Preparatory code: 
%------------------------------------------------------------
N_samples=length(Y);     N_ROC_points=N_samples+1;
N_patients=length(unique(patient_ID));
Positive_patients=unique( patient_ID(find(Y==1)) );
N_positive_patients=length( Positive_patients );
N_images=4*N_patients;
Tot_positives=length(find(Y==1));
Tot_negatives=N_samples-Tot_positives;

numOfFA=zeros(N_ROC_points,1); 
numOfD=numOfFA; numOfDetPatients=numOfFA;
numOfFA(1,1)=0; numOfD(1,1)=0; numOfDetPatients(1,1)=0;
Patients_detected_till_now=[];

%------------------------------------------------------------
% Sorting makes the algo work much faster
%------------------------------------------------------------
dat=sortrows([p, patient_ID, Y]);
dat=flipud(dat);

%get the indexes where each "uniquely scored" batch ends
[u_roc_scores, batch_idx]=unique(dat(:,1));
N_unique=size(u_roc_scores, 1);
temp=flipud([u_roc_scores, batch_idx]);
u_roc_scores=temp(:,1);
batch_idx=temp(:,2);
batch_idx=[0;batch_idx];
%------------------------------------------------------------
% Since the ROC curve is computed over patients and it depends on the
% ordering of the candidates with the same score, compute the average 
% over the permutation of the orderings
% Thanks to Peter Gehler for pointing this out.
%------------------------------------------------------------
maxIter=50;
numOfFA_total=zeros(N_ROC_points,maxIter);
numOfDetPatients_total=numOfFA_total;
numOfD_total=numOfFA_total;

for i=1:maxIter
    numOfFA=zeros(N_ROC_points,1); 
    numOfD=numOfFA; numOfDetPatients=numOfFA;
    numOfFA(1,1)=0; numOfD(1,1)=0; numOfDetPatients(1,1)=0;
    Patients_detected_till_now=[];

    %first, permute the ordering in every "batch"
    new_dat=zeros(size(dat));
    for k=2:N_unique+1
        data_batch=dat(batch_idx(k-1)+1:batch_idx(k),:);
        n=size(data_batch,1);
        idx=randperm(n);
        data_batch=data_batch(idx,:);
        new_dat(batch_idx(k-1)+1:batch_idx(k),:)=data_batch;
    end
    patient_ID=new_dat(:,2); Y=new_dat(:,3); 
    %------------------------------------------------------------
    % Now, just sweep through the data once to get all points in ROC
    %------------------------------------------------------------
    for j = 2 : N_ROC_points
        if (Y(j-1,1)==1)
		  % Another detection, but FP remains the same
            numOfD(j,1)=numOfD(j-1,1)+1;
            numOfFA(j,1)=numOfFA(j-1);
            Patients_detected_till_now=union(Patients_detected_till_now, patient_ID(j-1,1) );
            numOfDetPatients(j,1)=size(Patients_detected_till_now, 2);
        else
            % Patients_detected_till_now remains the same;
            numOfD(j,1)=numOfD(j-1,1);        
            numOfDetPatients(j,1)=numOfDetPatients(j-1,1);		  
            numOfFA(j,1)=numOfFA(j-1)+1;
        end
    end
    numOfD_total(:,i)=numOfD;
    numOfDetPatients_total(:,i)=numOfDetPatients;
    numOfFA_total(:,i)=numOfFA;
end
%------------------------------------------------------------
% And normalize the final results in the form needed
%------------------------------------------------------------
numOfFA=mean(numOfFA_total,2);
numOfD=mean(numOfD_total,2);
numOfDetPatients=mean(numOfDetPatients_total,2);


Pfa = numOfFA / Tot_negatives;
Pd = numOfD / Tot_positives;

FA_per_image = numOfFA / N_images;
Pd_patient_wise = numOfDetPatients / N_positive_patients;

index1=min(find(FA_per_image>=FA_low));
index2=max(find(FA_per_image<=FA_high));
if FA_per_image(end)<FA_high
	warning('AUC is incorrect because input parameter FA_high was higher than the maximum FA_per_image from the other parameters');
end
AUC=Pd_patient_wise(index1:(index2-1))' * (FA_per_image((index1+1):index2) - FA_per_image(index1:(index2-1)) );
return;

function n=rows(x);
n=size(x,1);
return;

%
%

